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Abstract We consider a new identity involving integrals and sums of Bessel 
functions. The identity provides new ways to evaluate integrals of products of 
two Bessel functions. The identity is remarkably simple and powerful since the 
summand and integrand are of exactly the same form and the sum converges 
to the integral relatively fast for most cases. A proof and numerical examples 
of the identity are discussed. 

1 Introduction 

The Newtonian kernel 

K(r,r') = —^ r,r'eM 3 , (1) 
\r — r | 

is ubiquitous in mathematical physics and is essential to an understanding of 
both gravitation and electrostatics [5]. It is central in classical mechanics [6 
but plays an equally important role in quantum mechanics [llj where it medi- 
ates the dominant two-particle interaction in electronic Schrodinger equations 
of atoms and molecules. 

Although the Newtonian kernel has many beautiful mathematical proper- 
ties, the fact that it is both singular and long-ranged is awkward and expen- 
sive from a computational point of view [8] and this has led to a great deal of 
research into effective methods for its treatment. Of the many schemes that 
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have been developed, Ewald partitioning [3] , multipole methods [7] and Fourier 
transform techniques |15j are particularly popular and have enabled the sim- 
ulation of large-scale particulate and continuous systems, even on relatively 
inexpensive computers. 

A recent alternative pfllSl fTBTll?] to these conventional techniques is to re- 
solve ([!]), a non-separable kernel, into a sum of products of one-body functions 

oo / OO / 

K(r,r') = J2 Yl YiMYimt/Wry) = J2 <Wr)<Wr') 

L—0 m— — l n,l= in— — I 

(2) 

where Yi m (r) is a real spherical harmonic P3J 14.30.2] of the angular part of 
three dimensional vector r, 

oo 

Ki(r,r) = Airj ^== dk, 



Ji (z) is a Bessel function of the first kind [2] 10.2.2], and r — \r\. The 
resolution ([2| is computationally useful because it decouples the coordinates 
r and r' and allows the two-body interaction integral 

E[p a ,Pb] = J J p a (r)K(r,r')p b (r')drdr', (3) 

between densities p a (r) and Pb(r) to be recast as 

oo oo I 

E[pa,Pb] = X! X! X! A nlmB n l m , (4) 
n=0 1=0 m=—l 

where A n i m is a one-body integral of the product of p a (r) and <j) n i m {r). If the 
one-body integrals can be evaluated efficiently and the sum converges rapidly, 
Q may offer a more efficient route to E[p a ,p b ] than ([3]). 
The key question is how best to obtain the Ki resolution 

oo 

Kt^r') = J2 K m{r)K n i(r'). 

n=0 

Previous attempts [4"ll51fT§] yielded complicated K n i whose practical utility is 
questionable but, recently, we have discovered the remarkable identity 

J u (at)J u (bt) ^ J v (ari)J v (bn) ._, 
dt=y £„ (5) 

where a,b £ [0,tt], v = 4, |, |, . . . and e n is defined by 

= n = (6) 

1, n > 1 1 j 
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This yields [12] the functions 



(/>nlm(r) = \ Ji+i/2{m)Y lm (r), 

V rn 1 

and these provide a resolution which is valid provided r < n. 

The aim of this Communication is to prove an extended version of the 
identity ^ and demonstrate its viability in approximating the integral of 
Bessel functions. 



2 Preliminaries 

The Bessel function of the first kind J v {£) is defined by [20] 3.1 (8)] 

(-1)" (z 
r(v + n+l)n\ V2 

n— 

It follows from that 
is an entire function of z and we have 

lim J v (z) I — I = ; r. 

Gauss' hypergeometric function is defined by PQ 2.1.2] 
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where (u) fc is the Pochhammer symbol (or rising factorial), given by 

(u) k = u{u + l)-- ■ (u + k- 1) . 

The series ^ converges absolutely for \z\ < 1 1, 2.1.1]. If Re (c — a — b) > 0, 
we have pQ 2.2.2] 

/ a ,6 \ _ r(c)r(c-a-b) 

2Fl { c r( c -a)r(c-by (9) 

Many special functions can be defined in terms of the hypergeometric func- 
tion. In particular, the Gegenbauer (or ultraspherical) polynomials C„ (t) are 
defined by [TO] 9.8.19] 



c<^) = ^(-T4 2A ;^) 



(10) 



with n € No and 

N = {0,1,...}. 
In the sequel, we will use the following Lemmas. 



Lemma 1 For k G No, we have 



x 2 i^ + H~ M + A > 2 ), N<i. 



Proof Using the formula [TJ 3.1.12] 



2a,2& x+l\ r(a + b+ §) T (§) f a,b 2 



2Fl [a + b + ^—)- r( a + \)r(b + \) 2Fl ( 



r(a)r(6) 



i 



2 



iii ( 10 ), we obtain 



,(m-2*) , , _ 2 2 "- 2 *~ 1 r fr, - 2fc + 1) r ( M - fc) f-k^-k 
J2k () ~ r(i-fc)r(2 M -4fc)(2fc)! 2 ^ | 



since r ^_ k ^ = for k = 0, 1, 

Applying Euler's transformation (TJ 2.2.7] 

2 f i c ;i =(i-i) 2^1 I c ;a; 



to (12 1, we get 



r o*-a*), , _ 2>- 2fc - 1 r( At -2fc + i)r( M -fc) 
2fc w r(|-*)r (2 M -4A)(2fe)i 

x(l-.x 2 )'^ +2fe 2 F 1 ^ + k >\- fl + k ;x 2 



and (111 follows since 

2 2 ^- 2fc - 1 r (// - 2fc + |) r (// - fc) = (fc + i- M ) fc 

J r(i-fc)r(2/^-4fc)(2fc)! ~~ M 
Lemma 2 Lei i/ie function h(x; a) be defined by 

h(x; a) = \ A t W (l - ^" 2fe) (f ), < x < a 

I 0, a < x < 7r 

where < a < 7T, 

_ (-l) fc (2fc)!r( A ,-2fc)2 2 ^ 2fc - 1 
fcW " a 2fc + 1 r(2 M -2fc) 

Re (ju) > 2fc — | and fc e N . 
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Then, h(x] a) can be represented by the Fourier cosine series 

h(x;a) = > e n7 j — cos (rax) , (14) 
«=o \2 an ) 

where e n was defined ira[6| 

Proof For Re(tr) > -§, a > and k € N , we have [201 3.32] 
l fe 

(i - ey-^ eg® cos (^) dt = ^U^±^l, h+2k (a) . (15) 



o 



(2jfe)!r(«r) (2a) 

Replacing cr = /i — 2fe and a = na in ( |15j ) , we obtain 

l 

— — = 2a / A£(a) (1 - rj C^fc ' (t) cos (nat) dt 



or 



j/f ("«) „2fc _ 2 

(|an)' J 



n ^fc _ I for. a \ CQS /\ ^ 




and the result follows. 

3 Main result 

The discontinuous integral 



I(a,b)=J JAat ] x JAbt) dt, 
o 

was investigated by Weber [21] , Sonine [T5] and Schafheitlin [T7] . They proved 
that [23 13.4 (2)] 



a 



X-v-li 



2 



b»r( u+t *- x+1 ) ( „+ m -a+i ^-g-A+i / 6 \a 



J (a,b) = ^V2fi 2 



2 A I> + l)r( A+ ^' v+1 



2 



(17) 
for 

Re(/i + + 1) > Re (A) > -1 (18) 
and < b < a. The corresponding expression for the case when < a < b, is 



obtained from (17) by interchanging a,b and also fj,, v. When a — b, we have 
[201 13.41 (2)] 

1 a) ~ 2 i r ( X + n-i' + l \ r f \ + v-fi+l \ r fi' + n + \ + V 
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provided that Re(^t + u + l) > Re (A) > 0. This result also follows from Gauss' 
summation formula ^ and ( fl7| . 

Theorem 1 If0<b<a<n, Re (fij > 2k - \ , Re (v) > - \, k e No, and 



J M (an) J„ (bn) / an 



2k 



then, 

Sk (a, b) = - r [ k+ % ^ 2 f 1 { ^ m - m + * . ^ y y (19) 



(2 \ ^ — o 
1 - fa J and integrating from to b, 

we get 
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i^ w «.)/(.-^r(i-5r"^ ; 



2X M -2fc-i 







(20) 

where we have used the integral representation ( 16 ). Setting x — bt and b = ua 
in pOl), we obtain 



o 

(21) 



Thus, we can re- write (21 1 as 



2b a» ™ a, ^ ^- 2k ~ l r (n - 2k + 1) r ( M - k) 



S k (a, b) = —Al (b) At (a) 

- 3 k r(| - fc) r (2/i - 4fc) (2fe)i 



7T 



x / (l-i 2 )^ f 2 Fr( ^ + H ^ + fc ;c^ 2 U, 
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or, using (13) and changing variables in the integral, 

{-l) k 2 2k ^ 



Sk (a, b) 



a ^r(u-k + \)r(v + \)r(\-k) 
i 

x Js-i(l-s) u -^ 2 F 1 (^ + k '\~ fl + k ;uj 2 S S jds. 



o 



Recalling the formula [TJ Th. 2.2.4] 



V 



valid for Re (c) > Re(d) > 0, x G C \ [1, oo) , we conclude that 



S k (a, b) = 



a ^r(a-k + \)r( v + \)r(\-k) 
r(i/ + i) 2 1 



2 

^+1 



But since [H 1.2.1] 



the result follows. 



The special case of Theorem [T] in which k = 0, was derived by Cooke in 
[2], as part of his work on Schlomilch series. 



Corollary 1 If < b < a < tt, Re > 2k - §, Re (v) > -|, k G N , th 



en 



J M (at) J v (U) 



-2 k 



dt = y^e n 



J M (an) J v (bn) 



n 



Proof The result follows immediately from (17) and (19), after taking A = 
u + v — 2k. Note that since for all k £ No 

Re(/x + v + 1) = Re (2k + 1 + A) > Re (1 + A) > Re (A) 



and 



Re (A) = Re (u + v - 2k) > -1, 



the conditions ( 18 ) are satisfied. 



Corollary 2 If < a,b < ir, and Re(/i) > 0, then 

oo 



J„ (an) J M (6n) _ f J M (at) J M (W) M _ j i (§)" , a < 6 



ra=0 



I 2/i V a/ 



a > 6 



Proof The result is a consequence of Corollary [T] and a special case of the 
integral §Tf\ (see [201 13.42 (1)]). 



s 




Fig. 2 Truncation error R.2o(o,,b) for a, b 6 [— 2tt,2tt] 
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2n 


Table 1 


i?jv(a, V) for 


a = 7r/2 and 


various b and JV 






N 


7r/4 


2tt/4 


3tt/4 


b 

An /A 


5ir/A 


6-k/A 


1 


3.15E-02 


1.81E-01 


3.15E-02 


-2.37E-02 


-4.12E-02 


-3.09E-2 


10 


-3.10E-03 


2.02E-02 


-2.40E-03 


1.45E-04 


1.74E-03 


-1.16E-2 


100 


-1.79E-06 


2.03E-03 


-4.81E-07 


-1.39E-07 


7.46E-07 


-1.17E-3 


1000 


1.81E-09 


2.03E-04 


4.86E-10 


-1.37E-10 


-7.46E-10 


-1.17E-4 


10000 


1.81E-12 


2.03E-05 


4.86E-13 


-1.37E-13 


-7.46E-13 


-1.17E-5 
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4 Numerical results 



The efficient and accurate computation of the 2 -Pi function in ([17]) can be 
numerically challenging [16] , and we note that Corollary 1 provides a new route 
to its evaluation. Because the numerator is bounded | J M (at) (bt) \ < 1 (2Ql 
13.42 (10)], the Weierstrass M-test shows that the series converges uniformly 
as long as /1 + v — 2k > 1 and fi, v > 0. 

To explore the rate of convergence of the sum in Corollary 2, we choose 
v = 5/2. The exact value of the integral is 



oc 

/ 



J5/2(at)J 5 / 2 (bt) _ l ^ 
t ~ 5 1 



min(|a|, \b\) 



-1 5/2 



(22) 



and truncation of the infinite series yields the finite sum 



J5/2 (an) Jc/2 (bn) 
T N (a, b) = V 5/ V ; 5/ y ' 



(23) 



and a truncation error 



R N {a, b) 



min(|q|, 
max(|a|, |6|) 



5/2 



T N (a,b). 



(24) 



The integral in (22) and sum in (23) are illustrated in Figure 1 and their 
difference (24) is shown in Figure 2 and Table 1. 

The truncation errors in Table 1 reveal that the rate of convergence of the 
Bessel sum is strongly dependent on the value of a and b. For b = 2n/A and 
b = 67r/4, the truncation error appears to decay as 1/N but, for the other 
values of 6, it appears to decay as 1/iV 3 . Although these empirical results are 
interesting, we have not been able to determine analytically the decay behavior 
as a function of a and b. 

The excellent agreement in Figure 1 and apparently flat plateau in Figure 
2 strongly suggest that Corollary 2 is true over the larger domain \a\ + \b\ < 2n. 
However, we have not yet managed to find a proof for this. 



5 Concluding remark 

We provide a rigorous proof that the identity (5) is valid on the square domain 
a, b E (0,7r). A numerical study indicates that the rate of convergence of the 
sum in the identity is sensitive to the values of a and b and further work to 
quantify this would be helpful. Generalization of the identity is possible and 
should be explored in the future work. 
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